<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Figure 6.24: Fitting a convex function to given data</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/convex_interpolation.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Figure 6.24: Fitting a convex function to given data</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.5.5</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Argyris Zymnis - 11/27/2005</span>
<span class="comment">%</span>
<span class="comment">% Here we find the convex function f that best fits</span>
<span class="comment">% some given data in the least squares sense.</span>
<span class="comment">% To do this we solve</span>
<span class="comment">%     minimize    ||yns - yhat||_2</span>
<span class="comment">%     subject to  yhat(j) &gt;= yhat(i) + g(i)*(u(j) - u(i)), for all i,j</span>

clear

<span class="comment">% Noise level in percent and random seed.</span>
rand(<span class="string">'state'</span>,29);
noiseint=.05;

<span class="comment">% Generate the data set</span>
u = [0:0.04:2]';
m=length(u);
y = 5*(u-1).^4 + .6*(u-1).^2 + 0.5*u;
v1=u&gt;=.2;
v2=u&lt;=.6;
v3=v1.*v2;
dipvec=((v3.*u-.4*ones(1,size(v3,2))).^(2)).*v3;
y=y+40*(dipvec-((.2))^2*v3);

<span class="comment">% add perturbation and plots the input data</span>
randf=noiseint*(rand(m,1)-.5);
yns=y+norm(y)*(randf);
figure
plot(u,yns,<span class="string">'o'</span>);

<span class="comment">% min. ||yns-yhat||_2</span>
<span class="comment">% s.t. yhat(j) &gt;= yhat(i) + g(i)*(u(j) - u(i)), for all i,j</span>
cvx_begin
    variables <span class="string">yhat(m)</span> <span class="string">g(m)</span>
    minimize(norm(yns-yhat))
    subject <span class="string">to</span>
        yhat*ones(1,m) &gt;= ones(m,1)*yhat' + (ones(m,1)*g').*(u*ones(1,m)-ones(m,1)*u');
cvx_end

nopts =1000;
t = linspace(0,2,nopts);
f = max(yhat(:,ones(1,nopts)) + <span class="keyword">...</span>
      g(:,ones(1,nopts)).*(t(ones(m,1),:)-u(:,ones(1,nopts))));
plot(u,yns,<span class="string">'o'</span>,t,f,<span class="string">'-'</span>);
axis <span class="string">off</span>
<span class="comment">%print -deps interpol_convex_function2.eps</span>
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Calling SDPT3: 2602 variables, 103 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 103
 dim. of socp   var  = 52,   num. of socp blk  =  1
 dim. of linear var  = 2550
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|1.1e+04|5.1e+01|3.4e+07|-2.463801e-08  0.000000e+00| 0:0:00| chol  1  1 
 1|0.023|0.867|1.1e+04|6.9e+00|3.2e+07| 7.051700e+01 -3.679632e+03| 0:0:00| chol  1  1 
 2|0.918|0.974|8.6e+02|2.7e-01|2.7e+06| 1.401794e+01 -2.943432e+03| 0:0:00| chol  1  1 
 3|0.910|0.987|7.8e+01|3.1e-02|2.4e+05| 1.029536e+01 -6.757387e+02| 0:0:00| chol  1  1 
 4|0.961|0.969|3.1e+00|9.2e-03|1.0e+04| 1.085902e+01 -4.145262e+02| 0:0:00| chol  1  1 
 5|0.965|0.985|1.1e-01|9.8e-04|6.9e+02| 1.092364e+01 -3.723299e+02| 0:0:00| chol  1  1 
 6|0.229|0.314|8.2e-02|7.0e-04|5.6e+02| 1.109744e+01 -3.294966e+02| 0:0:00| chol  1  1 
 7|0.219|0.312|6.4e-02|1.7e-02|4.9e+02| 1.114592e+01 -3.045759e+02| 0:0:00| chol  1  1 
 8|0.720|0.287|1.8e-02|2.5e-02|3.4e+02| 1.122282e+01 -2.740088e+02| 0:0:00| chol  1  1 
 9|0.739|0.306|4.7e-03|2.1e-02|2.6e+02| 1.125567e+01 -2.344734e+02| 0:0:00| chol  1  1 
10|0.217|0.246|3.6e-03|1.7e-02|2.3e+02| 1.127324e+01 -2.101866e+02| 0:0:00| chol  1  1 
11|0.187|0.582|3.0e-03|7.7e-03|1.9e+02| 1.132396e+01 -1.716350e+02| 0:0:00| chol  1  1 
12|0.343|1.000|1.9e-03|5.9e-04|1.4e+02| 1.135751e+01 -1.171362e+02| 0:0:00| chol  1  1 
13|0.911|1.000|1.7e-04|3.9e-04|8.5e+01| 1.130069e+01 -7.287643e+01| 0:0:00| chol  1  1 
14|0.574|1.000|7.4e-05|3.5e-05|5.7e+01| 1.105172e+01 -4.493234e+01| 0:0:00| chol  1  1 
15|0.722|1.000|2.1e-05|1.5e-05|3.7e+01| 1.070173e+01 -2.550582e+01| 0:0:00| chol  1  1 
16|0.507|1.000|1.0e-05|4.1e-06|2.9e+01| 9.626004e+00 -1.850222e+01| 0:0:00| chol  1  1 
17|0.798|0.987|2.1e-06|2.1e-06|1.8e+01| 8.940952e+00 -8.677965e+00| 0:0:00| chol  1  1 
18|0.306|1.000|1.4e-06|4.1e-07|2.5e+01| 3.963923e+00 -2.041665e+01| 0:0:00| chol  1  1 
19|0.820|0.965|2.6e-07|3.0e-07|8.7e+00| 3.284605e+00 -5.242031e+00| 0:0:00| chol  1  1 
20|0.777|0.882|5.7e-08|8.7e-08|5.3e+00| 5.203090e-01 -4.702256e+00| 0:0:00| chol  1  1 
21|0.517|0.735|2.8e-08|3.4e-08|3.5e+00|-1.037220e-01 -3.567379e+00| 0:0:00| chol  1  1 
22|0.798|1.000|5.6e-09|5.5e-09|1.6e+00|-1.131497e+00 -2.689229e+00| 0:0:00| chol  1  1 
23|0.608|0.919|2.2e-09|1.6e-09|9.1e-01|-1.575050e+00 -2.467360e+00| 0:0:00| chol  1  1 
24|0.844|1.000|3.4e-10|4.9e-10|3.3e-01|-2.030276e+00 -2.351939e+00| 0:0:00| chol  1  1 
25|0.748|0.880|8.6e-11|2.4e-10|1.2e-01|-2.186506e+00 -2.299985e+00| 0:0:00| chol  1  1 
26|0.842|1.000|1.4e-11|3.6e-10|3.5e-02|-2.247701e+00 -2.281961e+00| 0:0:00| chol  1  1 
27|0.896|1.000|1.4e-12|4.8e-10|4.1e-03|-2.270685e+00 -2.274766e+00| 0:0:00| chol  1  1 
28|0.952|0.961|2.2e-13|4.5e-10|2.0e-04|-2.274089e+00 -2.274289e+00| 0:0:00| chol  1  1 
29|0.922|0.986|3.2e-12|7.0e-10|4.0e-05|-2.274236e+00 -2.274276e+00| 0:0:00| chol  1  1 
30|0.658|0.991|1.2e-12|9.9e-10|1.5e-05|-2.274253e+00 -2.274268e+00| 0:0:00| chol  1  1 
31|0.920|1.000|6.0e-13|1.4e-09|2.0e-06|-2.274264e+00 -2.274266e+00| 0:0:00| chol  1  1 
32|0.807|1.000|3.5e-12|1.5e-09|4.4e-07|-2.274265e+00 -2.274266e+00| 0:0:00| chol  1  1 
33|0.890|1.000|8.0e-12|2.2e-09|8.1e-08|-2.274266e+00 -2.274266e+00| 0:0:00|
  stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 33
 primal objective value = -2.27426576e+00
 dual   objective value = -2.27426584e+00
 gap := trace(XZ)       = 8.06e-08
 relative gap           = 1.45e-08
 actual relative gap    = 1.44e-08
 rel. primal infeas     = 7.98e-12
 rel. dual   infeas     = 2.15e-09
 norm(X), norm(y), norm(Z) = 2.7e+00, 3.2e+08, 2.6e+09
 norm(A), norm(b), norm(C) = 8.4e+01, 2.0e+00, 1.3e+02
 Total CPU time (secs)  = 0.46  
 CPU time per iteration = 0.01  
 termination code       =  0
 DIMACS: 8.0e-12  0.0e+00  3.2e-08  0.0e+00  1.4e-08  1.5e-08
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +2.27427
 
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="convex_interpolation__01.png" alt=""> 
</div>
</div>
</body>
</html>